318 8.2 Molecular Simulation Methods
In other words, this probability distribution is a Gaussian curve with zero mean and vari
ance σv
B
i
k T m
2 =
/
. To determine the initial velocity, a pseudorandom probability from a uni
form distribution in the range 0–1 is generated for each atom and each spatial dimensional,
which is then sampled from the Maxwell–Boltzmann distribution to interpolate an equiva
lent speed (either positive or negative) parallel to x, y, and z for a given of set system tempera
ture T. The net momentum P is then calculated, for example, for the x component as
(8.3)
p
m v
x
i
n
i
x i
=
=∑
1
,
This is then used to subtract away a momentum offset from each initial atomic velocity so
that the net starting momentum of the whole system is zero (to minimize computational
errors due to large absolute velocity values evolving), so again for the x component
(8.4)
v
v
p
m
x i
x i
x
i
, ,
, ,
0
0
=
−
and similarly, for y and z components. A similar correction is usually applied to set the net
angular momentum to zero to prevent the simulated biomolecule from rotating, unless
boundary conditions (e.g., in the case of explicit solvation, see in the following text) make
this difficult to achieve in practice.
Then Newton’s second law (F = ma) is used to determine the acceleration a and ultimately
the velocity v of each atom, and thus where each will be, position vector of magnitude r, after
a given time (Figure 8.1), and to repeat this over a total simulation time that can be anything
from ~10 ps up to several hundred nanoseconds, usually imposed by a computational limit in
the absence of any coarse-graining to the simulations, in other words, to generate determin
istic trajectories of atomic positions as a function of time. This method involves numerical
integration over time, usually utilizing a basic algorithm called the “Verlet algorithm,” which
results in low round-off errors. The Verlet algorithm is developed from a Taylor expansion of
each atomic position:
(8.5)
r t
t
r t
v t
t
F t
m
t
F t
m
t
O
t
+
(
) = ( )+ ( )
+
( ) (
) +
( ) (
) +
(
)
(
)
∆
∆
∆
∆
∆
2
3
4
2
3!
Thus, the error is atomic position scales as ~O((Δt)4), which is small for the femtosecond time
increments normally employed, with the equivalent recursive relation for velocities having
errors that scale as ~O((Δt)2). The velocity form of the Verlet algorithm most commonly used
in MD can be summarized as
(8.6)
r t
t
r t
v t
t
a t
t
v t
t
v t
a t
+
(
) = ( )+ ( )
+ (
) ( )(
)
+
(
) = ( )+ (
) ( )
∆
∆
∆
∆
1 2
2
1 2
2
/
/
/
∆
∆
∆
∆
∆
∆
t
a t
t
m
U r t
t
v t
t
v t
t
a t
t
+
(
) = −(
)∇
=
(
)
(
)
+
(
) =
=
(
)+ (
) ( )
1
2
1 2
/
/
/
Related, though less popular, variants of this Verlet numerical integration include the
LeapFrog algorithm (equivalent to the Verlet integration but interpolates for velocities at
half time step intervals, but has a disadvantage in that velocities are not known at the same
time as the atomic positions, hence the name) and Beeman algorithm (generates identical
estimates for atomic positions as the Verlet method, but uses a slightly more accurate but
computationally more costly method to estimate velocities).